home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Amiga Plus 1995 #5 & #6
/
Amiga Plus CD - 1995 - No. 5 and 6.iso
/
pd
/
serien
/
purity
/
nr.49
/
astronomie
/
units
/
moolib.p
< prev
next >
Wrap
Text File
|
1995-06-24
|
20KB
|
473 lines
{---------------------------------------------------------------------------}
{ Unit MOOLIB: Mondbahn }
{---------------------------------------------------------------------------}
UNIT MooLib;
{$Projekt Astronomie}
INTERFACE;
FROM Astronomie USES MatLib,PnuLib,SphLib;
{---------------------------------------------------------------------------}
{ MINI_MOON: Mondkoordinaten geringer Genauigkeit (ca.5'/1') }
{ T : Zeit in jul.Jahrh. seit J2000 ( T=(JD-2451545)/36525 ) }
{ RA : Rektaszension (in h) }
{ DEC: Deklination (in Grad) }
{ (Aequinoktium des Datums) }
{---------------------------------------------------------------------------}
Procedure MINI_MOON(T:Real; Var RA,DEC:Real);
{---------------------------------------------------------------------------}
{ }
{ MOON: analytische Mondtheorie nach E.W.Brown (Improved Lunar Ephemeris) }
{ mit einer Genauigkeit von ca. 1" }
{ }
{ T: Zeit in julianischen Jahrhunderten seit J2000 (Ephemeridenzeit) }
{ (T=(JD-2451545.0)/36525.0) }
{ LAMBDA: geozentrische ekliptikale Laenge (Aequinoktium des Datums) }
{ BETA: geozentrische ekliptikale Breite (Aequinoktium des Datums) }
{ R: geozentrische Entfernung (in Erdradien) }
{ }
{---------------------------------------------------------------------------}
Procedure MOON(T:Real; Var LAMBDA,BETA,R: Real);
{---------------------------------------------------------------------------}
{ MOONEQU: aequatoriale Mondkoordinaten }
{ (Rektaszension RA und Deklination DEC in Grad, R in Erdradien) }
{ T in julian.Jahrhndt. seit J2000 ( T:= (JD - 2451545.0)/36525 ) }
{ Die Koord. beziehen sich auf das wahre Aequinoktium des Datums. }
{---------------------------------------------------------------------------}
Procedure MOONEQU(T:Real; Var RA,DEC,R: Real);
{---------------------------------------------------------------------------}
{ T_FIT_MOON: Berechnet die Tschebyscheff-Entwicklung der }
{ Koordinaten des Mondes (Reihen fuer RA,DEC und Radius). }
{ }
{ TA : Beginn des Entwicklungsintervalls (jul.Jahrh. seit J2000) }
{ TB : Ende des Entwicklungsintervalls ( TB < TA + 1 Monat ) }
{ N : Ordnung der Entwicklung }
{ RA_POLY,DE_POLY,R_POLY: Tschebyscheff Polynome fuer RA,DEC,R }
{---------------------------------------------------------------------------}
Procedure T_FIT_MOON(TA,TB:Real; N:Integer;
Var RA_POLY,DE_POLY,R_POLY:TPOLYNOM);
IMPLEMENTATION
{---------------------------------------------------------------------------}
{ MINI_MOON: Mondkoordinaten geringer Genauigkeit (ca.5'/1') }
{ T : Zeit in jul.Jahrh. seit J2000 ( T=(JD-2451545)/36525 ) }
{ RA : Rektaszension (in h) }
{ DEC: Deklination (in Grad) }
{ (Aequinoktium des Datums) }
{---------------------------------------------------------------------------}
Procedure MINI_MOON;
Const P2 = 6.283185307;
ARC = 206264.8062;
COSEPS = 0.91748;
SINEPS = 0.39778; { cos/sin(Ekliptikschiefe) }
Var L0,L,LS,F,D,H,S,N,DL,CB : Real;
L_MOON,B_MOON,V,W,X,Y,Z,RHO: Real;
Begin
{ mittlere Elemente der Mondbahn }
L0:= frac(0.606433+1336.855225*T); { mittl. Laenge des Mondes (in r) }
L :=P2*frac(0.374897+1325.552410*T); { mittl. Anomalie des Mondes }
LS:=P2*frac(0.993133+ 99.997361*T); { mittl. Anomalie Sonne }
D :=P2*frac(0.827361+1236.853086*T); { Diff. Laenge Mond-Sonne }
F :=P2*frac(0.259086+1342.227825*T); { Knotenabstand }
DL := +22640*sin(L) - 4586*sin(L-2*D) + 2370*sin(2*D) + 769*sin(2*L)
-668*sin(LS) - 412*sin(2*F) - 212*sin(2*L-2*D) - 206*sin(L+LS-2*D)
+192*sin(L+2*D) - 165*sin(LS-2*D) - 125*sin(D) - 110*sin(L+LS)
+148*sin(L-LS) - 55*sin(2*F-2*D);
S := F + (DL+412*sin(2*F)+541*sin(LS)) / ARC;
H := F-2*D;
N := -526*sin(H) + 44*sin(L+H) - 31*sin(-L+H) - 23*sin(LS+H)
+ 11*sin(-LS+H) -25*sin(-2*L+F) + 21*sin(-L+F);
L_MOON := P2 * frac ( L0 + DL/1296E3 ); { in rad }
B_MOON := ( 18520.0*sin(S) + N ) / ARC; { in rad }
{ aequatoriale Koordinaten }
CB := cos(B_MOON);
X := CB*cos(L_MOON);
V := CB*sin(L_MOON);
W := sin(B_MOON);
Y := COSEPS*V-SINEPS*W;
Z := SINEPS*V+COSEPS*W;
RHO:= sqrt(1.0-Z*Z);
DEC := (360.0/P2)*arctan(Z/RHO);
RA := ( 48.0/P2)*arctan(Y/(X+RHO));
if RA<0 then RA:=RA+24.0;
End;
{---------------------------------------------------------------------------}
{ }
{ MOON: analytische Mondtheorie nach E.W.Brown (Improved Lunar Ephemeris) }
{ mit einer Genauigkeit von ca. 1" }
{ }
{ T: Zeit in julianischen Jahrhunderten seit J2000 (Ephemeridenzeit) }
{ (T=(JD-2451545.0)/36525.0) }
{ LAMBDA: geozentrische ekliptikale Laenge (Aequinoktium des Datums) }
{ BETA: geozentrische ekliptikale Breite (Aequinoktium des Datums) }
{ R: geozentrische Entfernung (in Erdradien) }
{ }
{---------------------------------------------------------------------------}
Procedure MOON;
Const PI2 = 6.283185308;
ARC = 206264.81; { 3600*180/pi = Bogensekunden pro radian }
Var DGAM,FAC : Real;
DLAM,N,GAM1C,SINPI : Real;
L0, L, LS, F, D ,S : Real;
DL0,DL,DLS,DF,DD,DS: Real;
CO,SI: Array[-6..6,1..4] of Real;
{ berechne c=cos(a1+a2) und s=sin(a1+a2) aus den Additionstheo- }
{ remen fuer c1=cos(a1), s1=sin(a1), c2=cos(a2) und s2=sin(a2) }
Procedure ADDTHE(C1,S1,C2,S2:Real;Var C,S:Real);
Begin
C:=C1*C2-S1*S2;
S:=S1*C2+C1*S2;
End;
{ berechne sin(phi); phi in Einheiten von 1r=360 Grad }
Function SINUS(PHI:Real):Real;
Begin
SINUS:=sin(PI2*frac(PHI));
End;
{ berechne die langperiodischen Aenderungen der mittleren Elemente }
{ l,l',F,D und L0 sowie dgamma }
Procedure LONG_PERIODIC(T:Real; Var DL0,DL,DLS,DF,DD,DGAM:Real);
Var S1,S2,S3,S4,S5,S6,S7: Real;
Begin
S1:=SINUS(0.19833+0.05611*T);
S2:=SINUS(0.27869+0.04508*T);
S3:=SINUS(0.16827-0.36903*T);
S4:=SINUS(0.34734-5.37261*T);
S5:=SINUS(0.10498-5.37899*T);
S6:=SINUS(0.42681-0.41855*T);
S7:=SINUS(0.14943-5.37511*T);
DL0:= 0.84*S1+0.31*S2+14.27*S3+ 7.26*S4+ 0.28*S5+0.24*S6;
DL := 2.94*S1+0.31*S2+14.27*S3+ 9.34*S4+ 1.12*S5+0.83*S6;
DLS:=-6.40*S1 -1.89*S6;
DF := 0.21*S1+0.31*S2+14.27*S3-88.70*S4-15.30*S5+0.24*S6-1.86*S7;
DD := DL0-DLS;
DGAM := -3332E-9 * sinUS(0.59734-5.37261*T)
-539E-9 * sinUS(0.35498-5.37899*T)
-64E-9 * sinUS(0.39943-5.37511*T);
End;
{ INIT: berechne die mittleren Elemente und deren sin und cos }
{ l Anomalie des Mondes l' Anomalie der Sonne }
{ F Abstand des Mondes vom Knoten D Elongation des Mondes }
Procedure INIT;
Var I,J,MAX : Integer;
T2,ARG,FAC: Real;
Begin
T2:=T*T;
DLAM :=0; DS:=0; GAM1C:=0; sinPI:=3422.7000;
LONG_PERIODIC ( T, DL0,DL,DLS,DF,DD,DGAM );
L0 := PI2*frac(0.60643382+1336.85522467*T-0.00000313*T2) + DL0/ARC;
L := PI2*frac(0.37489701+1325.55240982*T+0.00002565*T2) + DL /ARC;
LS := PI2*frac(0.99312619+ 99.99735956*T-0.00000044*T2) + DLS/ARC;
F := PI2*frac(0.25909118+1342.22782980*T-0.00000892*T2) + DF /ARC;
D := PI2*frac(0.82736186+1236.85308708*T-0.00000397*T2) + DD /ARC;
for I := 1 to 4 do
Begin
Case I of
1: Begin ARG:=L; MAX:=4; FAC:=1.000002208; End;
2: Begin ARG:=LS; MAX:=3; FAC:=0.997504612-0.002495388*T; End;
3: Begin ARG:=F; MAX:=4; FAC:=1.000002708+139.978*DGAM; End;
4: Begin ARG:=D; MAX:=6; FAC:=1.0; End;
End;
CO[0,I]:=1.0;
CO[1,I]:=COS(ARG)*FAC;
SI[0,I]:=0.0;
SI[1,I]:=sin(ARG)*FAC;
for J := 2 to MAX do
ADDTHE(CO[J-1,I],SI[J-1,I],CO[1,I],SI[1,I],CO[J,I],SI[J,I]);
for J := 1 to MAX do
Begin
CO[-J,I]:=CO[J,I];
SI[-J,I]:=-SI[J,I];
End;
End;
End;
{ TERM berechne X=cos(p*arg1+q*arg2+r*arg3+s*arg4) und }
{ Y=sin(p*arg1+q*arg2+r*arg3+s*arg4) }
Procedure TERM(P,Q,R,S:Integer;Var X,Y:Real);
Var I: Array[1..4] of Integer; K: Integer;
Begin
I[1]:=P;
I[2]:=Q;
I[3]:=R;
I[4]:=S;
X :=1.0;
Y :=0.0;
for K:=1 to 4 do
if (I[K]<>0) Then ADDTHE(X,Y,CO[I[K],K],SI[I[K],K],X,Y);
End;
Procedure ADDSOL(COEFFL,COEFFS,COEFFG,COEFFP:Real;P,Q,R,S:Integer);
Var X,Y: Real;
Begin
TERM(P,Q,R,S,X,Y);
DLAM :=DLAM +COEFFL*Y;
DS :=DS +COEFFS*Y;
GAM1C:=GAM1C+COEFFG*X;
SINPI:=SINPI+COEFFP*X;
End;
Procedure SOLAR1;
Begin
ADDSOL( 13.902, 14.06,-0.001, 0.2607,0, 0, 0, 4);
ADDSOL( 0.403, -4.01,+0.394, 0.0023,0, 0, 0, 3);
ADDSOL( 2369.912, 2373.36,+0.601, 28.2333,0, 0, 0, 2);
ADDSOL( -125.154, -112.79,-0.725, -0.9781,0, 0, 0, 1);
ADDSOL( 1.979, 6.98,-0.445, 0.0433,1, 0, 0, 4);
ADDSOL( 191.953, 192.72,+0.029, 3.0861,1, 0, 0, 2);
ADDSOL( -8.466, -13.51,+0.455, -0.1093,1, 0, 0, 1);
ADDSOL(22639.500,22609.07,+0.079, 186.5398,1, 0, 0, 0);
ADDSOL( 18.609, 3.59,-0.094, 0.0118,1, 0, 0,-1);
ADDSOL(-4586.465,-4578.13,-0.077, 34.3117,1, 0, 0,-2);
ADDSOL( +3.215, 5.44,+0.192, -0.0386,1, 0, 0,-3);
ADDSOL( -38.428, -38.64,+0.001, 0.6008,1, 0, 0,-4);
ADDSOL( -0.393, -1.43,-0.092, 0.0086,1, 0, 0,-6);
ADDSOL( -0.289, -1.59,+0.123, -0.0053,0, 1, 0, 4);
ADDSOL( -24.420, -25.10,+0.040, -0.3000,0, 1, 0, 2);
ADDSOL( 18.023, 17.93,+0.007, 0.1494,0, 1, 0, 1);
ADDSOL( -668.146, -126.98,-1.302, -0.3997,0, 1, 0, 0);
ADDSOL( 0.560, 0.32,-0.001, -0.0037,0, 1, 0,-1);
ADDSOL( -165.145, -165.06,+0.054, 1.9178,0, 1, 0,-2);
ADDSOL( -1.877, -6.46,-0.416, 0.0339,0, 1, 0,-4);
ADDSOL( 0.213, 1.02,-0.074, 0.0054,2, 0, 0, 4);
ADDSOL( 14.387, 14.78,-0.017, 0.2833,2, 0, 0, 2);
ADDSOL( -0.586, -1.20,+0.054, -0.0100,2, 0, 0, 1);
ADDSOL( 769.016, 767.96,+0.107, 10.1657,2, 0, 0, 0);
ADDSOL( +1.750, 2.01,-0.018, 0.0155,2, 0, 0,-1);
ADDSOL( -211.656, -152.53,+5.679, -0.3039,2, 0, 0,-2);
ADDSOL( +1.225, 0.91,-0.030, -0.0088,2, 0, 0,-3);
ADDSOL( -30.773, -34.07,-0.308, 0.3722,2, 0, 0,-4);
ADDSOL( -0.570, -1.40,-0.074, 0.0109,2, 0, 0,-6);
ADDSOL( -2.921, -11.75,+0.787, -0.0484,1, 1, 0, 2);
ADDSOL( +1.267, 1.52,-0.022, 0.0164,1, 1, 0, 1);
ADDSOL( -109.673, -115.18,+0.461, -0.9490,1, 1, 0, 0);
ADDSOL( -205.962, -182.36,+2.056, +1.4437,1, 1, 0,-2);
ADDSOL( 0.233, 0.36, 0.012, -0.0025,1, 1, 0,-3);
ADDSOL( -4.391, -9.66,-0.471, 0.0673,1, 1, 0,-4);
End;
Procedure SOLAR2;
Begin
ADDSOL( 0.283, 1.53,-0.111, +0.0060,1,-1, 0,+4);
ADDSOL( 14.577, 31.70,-1.540, +0.2302,1,-1, 0, 2);
ADDSOL( 147.687, 138.76,+0.679, +1.1528,1,-1, 0, 0);
ADDSOL( -1.089, 0.55,+0.021, 0.0 ,1,-1, 0,-1);
ADDSOL( 28.475, 23.59,-0.443, -0.2257,1,-1, 0,-2);
ADDSOL( -0.276, -0.38,-0.006, -0.0036,1,-1, 0,-3);
ADDSOL( 0.636, 2.27,+0.146, -0.0102,1,-1, 0,-4);
ADDSOL( -0.189, -1.68,+0.131, -0.0028,0, 2, 0, 2);
ADDSOL( -7.486, -0.66,-0.037, -0.0086,0, 2, 0, 0);
ADDSOL( -8.096, -16.35,-0.740, 0.0918,0, 2, 0,-2);
ADDSOL( -5.741, -0.04, 0.0 , -0.0009,0, 0, 2, 2);
ADDSOL( 0.255, 0.0 , 0.0 , 0.0 ,0, 0, 2, 1);
ADDSOL( -411.608, -0.20, 0.0 , -0.0124,0, 0, 2, 0);
ADDSOL( 0.584, 0.84, 0.0 , +0.0071,0, 0, 2,-1);
ADDSOL( -55.173, -52.14, 0.0 , -0.1052,0, 0, 2,-2);
ADDSOL( 0.254, 0.25, 0.0 , -0.0017,0, 0, 2,-3);
ADDSOL( +0.025, -1.67, 0.0 , +0.0031,0, 0, 2,-4);
ADDSOL( 1.060, 2.96,-0.166, 0.0243,3, 0, 0,+2);
ADDSOL( 36.124, 50.64,-1.300, 0.6215,3, 0, 0, 0);
ADDSOL( -13.193, -16.40,+0.258, -0.1187,3, 0, 0,-2);
ADDSOL( -1.187, -0.74,+0.042, 0.0074,3, 0, 0,-4);
ADDSOL( -0.293, -0.31,-0.002, 0.0046,3, 0, 0,-6);
ADDSOL( -0.290, -1.45,+0.116, -0.0051,2, 1, 0, 2);
ADDSOL( -7.649, -10.56,+0.259, -0.1038,2, 1, 0, 0);
ADDSOL( -8.627, -7.59,+0.078, -0.0192,2, 1, 0,-2);
ADDSOL( -2.740, -2.54,+0.022, 0.0324,2, 1, 0,-4);
ADDSOL( 1.181, 3.32,-0.212, 0.0213,2,-1, 0,+2);
ADDSOL( 9.703, 11.67,-0.151, 0.1268,2,-1, 0, 0);
ADDSOL( -0.352, -0.37,+0.001, -0.0028,2,-1, 0,-1);
ADDSOL( -2.494, -1.17,-0.003, -0.0017,2,-1, 0,-2);
ADDSOL( 0.360, 0.20,-0.012, -0.0043,2,-1, 0,-4);
ADDSOL( -1.167, -1.25,+0.008, -0.0106,1, 2, 0, 0);
ADDSOL( -7.412, -6.12,+0.117, 0.0484,1, 2, 0,-2);
ADDSOL( -0.311, -0.65,-0.032, 0.0044,1, 2, 0,-4);
ADDSOL( +0.757, 1.82,-0.105, 0.0112,1,-2, 0, 2);
ADDSOL( +2.580, 2.32,+0.027, 0.0196,1,-2, 0, 0);
ADDSOL( +2.533, 2.40,-0.014, -0.0212,1,-2, 0,-2);
ADDSOL( -0.344, -0.57,-0.025, +0.0036,0, 3, 0,-2);
ADDSOL( -0.992, -0.02, 0.0 , 0.0 ,1, 0, 2, 2);
ADDSOL( -45.099, -0.02, 0.0 , -0.0010,1, 0, 2, 0);
ADDSOL( -0.179, -9.52, 0.0 , -0.0833,1, 0, 2,-2);
ADDSOL( -0.301, -0.33, 0.0 , 0.0014,1, 0, 2,-4);
ADDSOL( -6.382, -3.37, 0.0 , -0.0481,1, 0,-2, 2);
ADDSOL( 39.528, 85.13, 0.0 , -0.7136,1, 0,-2, 0);
ADDSOL( 9.366, 0.71, 0.0 , -0.0112,1, 0,-2,-2);
ADDSOL( 0.202, 0.02, 0.0 , 0.0 ,1, 0,-2,-4);
End;
Procedure SOLAR3;
Begin
ADDSOL( 0.415, 0.10, 0.0 , 0.0013,0, 1, 2, 0);
ADDSOL( -2.152, -2.26, 0.0 , -0.0066,0, 1, 2,-2);
ADDSOL( -1.440, -1.30, 0.0 , +0.0014,0, 1,-2, 2);
ADDSOL( 0.384, -0.04, 0.0 , 0.0 ,0, 1,-2,-2);
ADDSOL( +1.938, +3.60,-0.145, +0.0401,4, 0, 0, 0);
ADDSOL( -0.952, -1.58,+0.052, -0.0130,4, 0, 0,-2);
ADDSOL( -0.551, -0.94,+0.032, -0.0097,3, 1, 0, 0);
ADDSOL( -0.482, -0.57,+0.005, -0.0045,3, 1, 0,-2);
ADDSOL( 0.681, 0.96,-0.026, 0.0115,3,-1, 0, 0);
ADDSOL( -0.297, -0.27, 0.002, -0.0009,2, 2, 0,-2);
ADDSOL( 0.254, +0.21,-0.003, 0.0 ,2,-2, 0,-2);
ADDSOL( -0.250, -0.22, 0.004, 0.0014,1, 3, 0,-2);
ADDSOL( -3.996, 0.0 , 0.0 , +0.0004,2, 0, 2, 0);
ADDSOL( 0.557, -0.75, 0.0 , -0.0090,2, 0, 2,-2);
ADDSOL( -0.459, -0.38, 0.0 , -0.0053,2, 0,-2, 2);
ADDSOL( -1.298, 0.74, 0.0 , +0.0004,2, 0,-2, 0);
ADDSOL( 0.538, 1.14, 0.0 , -0.0141,2, 0,-2,-2);
ADDSOL( 0.263, 0.02, 0.0 , 0.0 ,1, 1, 2, 0);
ADDSOL( 0.426, +0.07, 0.0 , -0.0006,1, 1,-2,-2);
ADDSOL( -0.304, +0.03, 0.0 , +0.0003,1,-1, 2, 0);
ADDSOL( -0.372, -0.19, 0.0 , -0.0027,1,-1,-2, 2);
ADDSOL( +0.418, 0.0 , 0.0 , 0.0 ,0, 0, 4, 0);
ADDSOL( -0.330, -0.04, 0.0 , 0.0 ,3, 0, 2, 0);
End;
{ Stoerungsanteil N der ekliptikalen Breite }
Procedure SOLARN(Var N:Real);
Var X,Y: Real;
Procedure ADDN(COEFFN:Real;P,Q,R,S:Integer);
Begin
TERM(P,Q,R,S,X,Y);
N:=N+COEFFN*Y
End;
Begin
N := 0.0;
ADDN(-526.069, 0, 0,1,-2);
ADDN( -3.352, 0, 0,1,-4);
ADDN( +44.297,+1, 0,1,-2);
ADDN( -6.000,+1, 0,1,-4);
ADDN( +20.599,-1, 0,1, 0);
ADDN( -30.598,-1, 0,1,-2);
ADDN( -24.649,-2, 0,1, 0);
ADDN( -2.000,-2, 0,1,-2);
ADDN( -22.571, 0,+1,1,-2);
ADDN( +10.985, 0,-1,1,-2);
End;
{ Stoerungen der ekliptikalen Laenge durch Venus und Jupiter }
Procedure PLANETARY(Var DLAM:Real);
Begin
DLAM := DLAM
+0.82*sinUS(0.7736 -62.5512*T)+0.31*sinUS(0.0466 -125.1025*T)
+0.35*sinUS(0.5785 -25.1042*T)+0.66*sinUS(0.4591+1335.8075*T)
+0.64*sinUS(0.3130 -91.5680*T)+1.14*sinUS(0.1480+1331.2898*T)
+0.21*sinUS(0.5918+1056.5859*T)+0.44*sinUS(0.5784+1322.8595*T)
+0.24*sinUS(0.2275 -5.7374*T)+0.28*sinUS(0.2965 +2.6929*T)
+0.33*sinUS(0.3132 +6.3368*T);
End;
Begin
INIT;
SOLAR1;
SOLAR2;
SOLAR3;
SOLARN(N);
PLANETARY(DLAM);
LAMBDA := 360.0*frac( (L0+DLAM/ARC) / PI2 );
S := F + DS/ARC;
FAC := 1.000002708+139.978*DGAM;
BETA := ( FAC*(18518.511+1.189+GAM1C)*sin(S)-6.24*sin(3*S)+N ) / 3600.0;
SINPI := SINPI * 0.999953253;
R := ARC / SINPI;
End;
{---------------------------------------------------------------------------}
{ MOONEQU: aequatoriale Mondkoordinaten }
{ (Rektaszension RA und Deklination DEC in Grad, R in Erdradien) }
{ T in julian.Jahrhndt. seit J2000 ( T:= (JD - 2451545.0)/36525 ) }
{ Die Koord. beziehen sich auf das wahre Aequinoktium des Datums. }
{---------------------------------------------------------------------------}
Procedure MOONEQU;
Var L,B,X,Y,Z: Real;
Begin
MOON(T,L,B,R); { ekliptikale Moondkoordinaten }
CART(R,B,L,X,Y,Z); { (mittleres Aequinoktium des Datums) }
ECLEQU(T,X,Y,Z); { Umwandlung in aequatoriale Koordinaten }
NUTEQU(T,X,Y,Z); { Nutation }
POLAR(X,Y,Z,R,DEC,RA);
End;
{---------------------------------------------------------------------------}
{ T_FIT_MOON: Berechnet die Tschebyscheff-Entwicklung der }
{ Koordinaten des Mondes (Reihen fuer RA,DEC und Radius). }
{ }
{ TA : Beginn des Entwicklungsintervalls (jul.Jahrh. seit J2000) }
{ TB : Ende des Entwicklungsintervalls ( TB < TA + 1 Monat ) }
{ N : Ordnung der Entwicklung }
{ RA_POLY,DE_POLY,R_POLY: Tschebyscheff Polynome fuer RA,DEC,R }
{---------------------------------------------------------------------------}
Procedure T_FIT_MOON;
Const NDIM = 27;
Var I,J,K : Integer;
FAC,BPA,BMA,PHI: Real;
T,H,RA,DE,R : Array[0..NDIM] of Real;
Begin
if (NDIM<2*MAX_TP_DEG+1) Then WriteLn(' NDIM zu klein in T_FIT_MOON');
if (N>MAX_TP_DEG) Then WriteLn(' N zu gross in T_FIT_MOON');
RA_POLY.M := N;
DE_POLY.M := N;
R_POLY.M := N;
RA_POLY.A := TA;
DE_POLY.A := TA;
R_POLY.A := TA;
RA_POLY.B := TB;
DE_POLY.B := TB;
R_POLY.B := TB;
BMA := (TB-TA)/2.0;
BPA := (TB+TA)/2.0;
FAC := 2.0/(N+1);
PHI := PI/(2*N+2); { h(k)=cos(pi*k/N/2) }
H[0] := 1.0;
H[1] := COS(PHI);
for I:=2 to (2*N+1) do
H[I] := 2*H[1]*H[I-1]-H[I-2];
for K:=1 to N+1 DO
T[K] := H[2*K-1]*BMA+BPA; { Stuetzstellen }
for K:=1 to N+1 DO
MOONEQU(T[K],RA[K],DE[K],R[K]);
for K := 2 TO N+1 DO { RA stetig machen in }
if (RA[K-1]<RA[K]) Then RA[K]:=RA[K]-360.0; { [-360,+360] !!!! }
for J := 0 to N do { Tscheb.-Koeffizienten }
Begin { C(j) berechnen }
PHI := PI*J/(2*N+2);
H[1] := COS(PHI);
for I:=2 to (2*N+1) do
H[I] := 2*H[1]*H[I-1]-H[I-2];
RA_POLY.C[J] := 0.0;
DE_POLY.C[J] := 0.0;
R_POLY.C[J] := 0.0;
for K:=1 to N+1 do
Begin
RA_POLY.C[J] := RA_POLY.C[J] + H[2*K-1]*RA[K];
DE_POLY.C[J] := DE_POLY.C[J] + H[2*K-1]*DE[K];
R_POLY.C[J] := R_POLY.C[J] + H[2*K-1]*R[K];
End;
RA_POLY.C[J]:=RA_POLY.C[J]*FAC;
DE_POLY.C[J]:=DE_POLY.C[J]*FAC;
R_POLY.C[J] :=R_POLY.C[J]*FAC;
End;
End;
BEGIN
END.